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We compute the effective dispersion and vibrational density of states (DOS) of two-dimensional 
sub-regions of three dimensional face centered cubic (FCC) crystals using both a direct projection- 
inversion technique and a Monte Carlo simulation based on a common underlying Hamiltonian. 
We study both a (111) and (100) plane. We show that for any given direction of wavevector, both 
(111) and (100) show an anomalous ^ q regime at low q where cj^ is the energy associated 
with the given mode and q is its wavenumber. The uo^ ^ q scaling should be expected to give 
rise to an anomalous DOS, Do;, at low u: ^ cj^ rather than the conventional Debye result: 
Duj ~ Ct;^. The DOS for (100) looks to be consistent with D^o ^ cj"^, while (111) shows something 
closer to the conventional Debye result at the smallest frequencies. In addition to the direct 
projection-inversion calculation, we perform Monte Carlo simulations to study the effects of finite 
sampling statistics. We show that finite sampling artifacts act as an effective disorder and bias 
Duj^ giving a behavior closer to D^j cj^ than Duj ^ (jO^. These results should have an important 
impact on the interpretation of recent studies of colloidal solids where the two-point displacement 
correlations can be obtained directly in real-space via microscopy. 



PACS codes: 82.70.Dd, 63.20.dd, 63.22.-m. 

I. INTRODUCTION 

Colloidal suspensions of spherical particles exhibit 
similar behavior to conventional condensed matter sys- 
tems [1 . Direct observation of colloidal particle tra- 
jectories using optical microscopy has been employed to 
study condensed matter phenomena like: crystal nucle- 
ation [2 , impurity frustrated crystallization [3 , glass 
transitions [4], melting |5 etc. When quenched into a 
solid-like state, particles fluctuate about their equilib- 
rium positions [BHS]. At long wavelength these fluctua- 
tions should be governed by some effective elasticity and 
one should be able to extract the corresponding moduli 
from observed displacement fluctuations. 

There are several ways to extract moduli from exper- 
imentally observed displacement fluctuations. Zahn et. 
al. [To] obtained bulk elastic constants by studying strain 
fluctuations in sub-regions of a 2D hexagonal colloidal 
crystal of increasing size and extrapolated to the infinite 
system size limit, von Griinberg et al. [11 performed 
a planewave decomposition of particle displacements to 
compute the dispersion and used it to extract long wave- 
length elastic constants. However these methods only 
work in spatially homogeneous systems. Much like light 
scattering methods [121 US] i they average out any spatial 
disorder. In disordered systems, like structural glasses or 
geometrically ordered systems with heterogeneous inter- 
actions [14| [15] , the spatial heterogeneity and its impact 
on the fluctuation spectrum is of central interest, so other 
methods must be employed. 

To study disordered systems, a third method has been 
employed by several groups [7'-^, T6^ . This approach ex- 
ploits a connection between the linear elastic response 
function and displacement covariance matrix for a sys- 
tem in thermal equilibrium: Giaj(3 = {uiaUj(3) /k^T [17\ 



where k}y is Boltzmann's constant and T is the temper- 
ature. Here Uia is the displacement of the i-th particle 
in the a-th direction and the angle brackets mean ther- 
mal average. In glassy colloids it has been used to study: 
the DOS [SldH], the so-called Boson peak, an excess in 
the DOS when compared to the Debye theory, [8] and 
the connection between structural relaxation and low fre- 
quency normal modes [191 [20]. For colloidal crystals, the 
method has been used to characterize the spectrum of 
normal modes and DOS [3 [21]. However this wealth of 
information comes at a price and requires good statistical 
estimates for all entries of Giajp- This raises a number of 
problems related to statistical convergence [16 that can 
lead to qualitative misinterpretation of the data from fi- 
nite samples, as we will show below. 

There is another serious difficulty to which this method 
is susceptible when applied to a three-dimensional (3D) 
colloidal system. To assemble Giaj(3^ one needs many 
uncorrelated snapshots of the system. In a typical mi- 
croscopy experiment, one has access to a single plane of 
the material and can not make simultaneous observations 
of particles above or below. However the observed in- 
plane fluctuations in the patch are mediated by degrees 
of freedom off the imaging plane that are not observable. 

Thus one is faced with two problems. First, it is 
not obvious what effect the embedding medium has on 
the Green's function assembled from measurements made 
only in a lower dimensional patch. This difficulty arises 
even in the case of a perfect, homogeneous crystal. Sec- 
ond, it is not clear how finite statistics bias the spectrum, 
especially in the presence of the restricted window of ob- 
servation. 

Recently Maggs and coworkers have made good 
progress on the first question [22H25] . They have shown, 
most dramatically, that the expected energy of plane 
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wave fluctuations in the restricted 2D patch is hnear 
in wavenumber rather than quadratic as is the case for 
systems governed by a Laplace- hke operator [26]. How- 
ever, their approach starts with a long wavelength cu- 
bic elasticity tensor and does not treat short wavelength 
fluctuations near the Brillouin Zone (BZ) boundary in 
a way that takes into account the structure of the un- 
derlying lattice. This makes it difficult to separate arti- 
facts in the DOS near van Hove peaks into those arising 
from the observations being restricted to a 2D patch and 
those arising from flnite sampling. Furthermore, in the 
continuum elasticity approach, the anisotropic compo- 
nent of the elastic modulus tensor must be set by hand 
rather than emerging naturally from the underlying lat- 
tice structure and particle interactions. 

In this work, for a perfectly homogeneous face centered 
cubic (FCC) crystal with harmonic, nearest neighbor in- 
teractions, we flrst examine the dispersion relation and 
normal mode structure of the Green's function assem- 
bled from observing only in-plane fluctuations of a crys- 
tal patch. Then we investigate the artifacts induced by 
flnite sampling statistics. 

In contrast to the approach of Maggs and coworkers, 
we work in a fully atomistic framework. This allows us to 
produce realistic results for the entire Brilloiun zone and 
calculate the exact dispersion relation, frequency spec- 
trum and normal modes for flnite 2D patches cut from 
(100) and (111) FCC crystal planes. Our dispersion re- 
sults agree in the long wavelength limit with those in 
refs [22 , 23 . However, in the rest of the Brillouin zone we 
flnd signiflcant qualitative differences between the longi- 
tudinal and transverse branches of the dispersion. More- 
over, in the low frequency portion of the [111] patch DOS, 
we see pronounced deviations from the scaling expected 
from the long wavelength dispersion relation. 

At the same time we perform a Monte- Carlo simula- 
tion using precisely the same Hamiltonian to make a di- 
rect, quantitative comparison between the true projected 
spectrum and those spectra inferred from flnite statistical 
samples. We flnd that flnite statistics induces important 
artifacts in the DOS: it smears out and shifts the van 
Hove singularities to lower frequency parts of the spec- 
trum in a manner similar to disorder [15 . This alters the 
inferred low frequency scaling behavior in a qualitative 
way. 



II. ELASTICITY OF A TWO DIMENSIONAL 
SUBSYSTEM 

We analytically compute the elastic Green's function 
of a 2-dimensional patch of atoms embedded in homoge- 
neous, face-centered-cubic (FCC) crystal with cubic pe- 
riodic boundaries. The domain and range of the Green's 
function of a patch, denoted by ^, of M atoms is a M?^ 
subspace of the conflguration space of the FCC crys- 
tal (having N atoms). The embedding FCC crystal is 
governed by a harmonic Hamiltonian having only, pair- 



wise interactions, connecting nearest neighbors with un- 
stressed springs. We denote the nearest neighbor spacing 
by a and measure all lengths as units of a. For this system 
the bond stiffness tensor is: Biaj(3 = KrijcxTij(3 [26 . Here 
K is the bond strength; Latin indices i, j label atoms, 
while Greek indices a, {3 label Cartesian axes. Vij is the 
unit vector pointing from atom i to one of its nearest 
neighbors, atom j. The non-zero entries of the Hessian 
matrix, H, are given by: Hiajf3 = Biajf3 when i ^ and 
Hiaj(3 = - ^jj^i Biaj(3 for i = j. 

To obain the Green's function, ^, of a P x P atom 
patch, we start by computing the full Green's function 
G = EI~^, of the embedding cubic, periodic crystal with 
C 4- atom cubic unit cells along each edge. The com- 
putation of G nominally amounts to a computationally 
expensive, inversion of the operator H. However, due to 
periodicity, H is a block-circulant matrix, and we use a 
Fourier space approach (detailed in appendix [A| which 
permits inversion in linear computer time and space. In 
practice we have used this method to invert Hessians of 
roughly 8 million particle systems in a few tens of minutes 
of computer time on a single workstation. Projecting G 
onto the observed degrees of freedom in the patch, we 
immediately obtain Q. 

We next discuss various properties (dispersion, mode 
structure and DOS) of Q for patches cut from the (111) 
and (100) planes of an FCC crystal. The embedding 
crystal has C cubic unit cell cells along each edge, while 
the patch has P atoms along each edge of its boundary. 
Figure [l] has an illustration of our patches. We take the 
(111) patch to be a triangle lattice with rhombic bound- 
aries and the (100) patch to be a square lattice with 
square bounds. 

A. Pseudo-dispersion 

For each vector q that lies in the flrst BZ of the patch's 
reciprocal lattice, planewave excitations are given by: 
i^iaLiT){q) = PaLiT)e''^''^'; Pa IS the Component of the 
polarization vector along the a-axis. Here the subscripts 
L and T label longitudinal and transverse polarizations 
respectively, while fi labels the equilibrium position of 
atom i. Because of boundary effects, plane waves are not, 

FIG. 1: Left: FCC crystal with C = 3 cubic unit cells along 
each edge. Center: Rhombic (111) FCC patch with P = 4 
atoms along each edge. Right: Square (100) FCC patch with 
P = 4 atoms along each edge. 
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strictly speaking, eigenmodes of Q even for high symme- 
try directions. Nevertheless, we compute the longitudinal 
(transverse) pseudo-dispersion oj'^{q) as: 



V^i(T)(?) 



'^iaL{T)Giaj(3'^ji3L{T) 

iaj(3 



(1) 



In all subsequent results, we report energies in units of 

An approximation for scalar elasticity governed by the 
Laplace operator in an infinite continuum [25 as well as 
later work for continuum elasticity governed by a long 
wavelength cubic elasticity tensor [22l |23] predict the 
anomalous non-linear dispersion of oo'^ ^ q (at low q) 
for the patch Green's function Q. Our approach differs 
from this previous work in that the underlying lattice and 
interactions dictate the symmetry of the long wavelength 
elastic modulus tensor and determine the behavior of the 
system at the BZ boundary - in particular the magnitude 
and location of the van Hovepeaks - in a realistic way 
for an FCC crystal. In figure [2] we show 00'^ /q^ for various 
patches cut from the (111) plane and for longitudinal and 
transverse planewaves traveling along an edge of the BZ. 
Note that a patch with P = 32 atoms along an edge cut 
from a crystal of size C = 32 is roughly half the size of 
the largest (111) plane that fits in the embedding FCC 
crystal. Figure |2] shows that patches of size P = 32 are 
already quite insensitive to the size of the full embed- 
ding crystal once the embedding crystal is just slightly 
larger than the patch itself. We can also see that the 
ratio of the longitudinal to transverse branch is not 3 as 
it would be for a true 2D triangular lattice with central 
force interactions. 

Figure |3] shows the complete dispersion scaled by the 
wavenumber for wavevectors q that lie in the first BZ of 
a (111) patch (P = 64, C = 64). We see that at low 
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FIG. 3: Computed dispersion scaled by wavenumber, 
/Ka^q, for a (111) patch (P = 64, C = 64) for waves in the 
first Brillouin zone. Top: Longitudinal. Bottom: Transverse. 
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FIG. 2: Computed longitudinal and transverse dispersion 
curves scaled by wavenumber for planewaves traveling along 
the next-nearest neighbor direction and lying in the first Bril- 
louin zone for (111) various patches. 



q^ there is little variation in the scaled dispersion along 
any particular radial direction in the BZ. The angular de- 
pendence of the dispersion has hexagonal symmetry, as 
it must. Moreover at low wavenumbers the dispersion is 
isotropic - particularly in the longitudinal case. A similar 
calculation indicates that the dispersion of a 2D triangu- 
lar lattice has a similar angular dependence as that of the 
(111) FCC patch. Note that a slight breaking of hexag- 
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FIG. 4: Computed longitudinal and transverse dispersion 
curves scaled by wavenumber for planewaves lying in the first 
Brillouin zone of various (100) patches. Top: For planewaves 
travelling in the nearest neighbor direction ([010] direction). 
Bottom: For planewaves travelling in the next-nearest neigh- 
bor direction ([Oil] direction). 



onal symmetry is apparent in figure [3] near the center of 
the BZ. This is an artifact of cutting a rhombic (111) 
patch from a cubic FCC crystal: the edges of the rhom- 
bus are not crystallographically equivalent to its shorter 
diagonal. 

The dispersion for square (100) patches essentially fol- 
lows the same non-linear scaling cu'^ ^ q for low but 
with significant differences compared to (111). Figure |4] 
shows the dispersion relation for various (100) patches 
and for planewaves traveling along the [010] (top) and 
[Oil] (bottom) directions. We see that longitudinal and 
transverse planewaves traveling along a diagonal of the 
square ([Oil] direction) are nearly degenerate (as in the 
Hard-Sphere Monte Carlo simulations of reference [23 ), 
while the planewaves along the the [010] direction are not. 
This suggests that the (100) patch dispersion is strongly 
anisotropic. Figure |5] shows the dispersion scaled by the 
wavenumber for wavevectors q that lie in the first BZ of 
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FIG. 5: Computed dispersion scaled by wavenumber: 
uj'^/Ka^q, for a (100) patch (P = 64, C = 64) for waves in the 
first Brillouin zone. Top: Longitudinal. Bottom: Transverse. 



a (100) patch (P = 64, C = 64). The scaled dispersion 
shows strong anisotropy that is particularly pronounced 
in the transverse branch. For the transverse case, uo'^ /q 
also shows radial variations at low q along some, but not 
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all, directions indicating that the low-g scaling is not as 
robust in every direction, with the [010] direction giving 
particularly poor scaling at small q. 



B. Density of States 

The Debye prediction for the low frequency scaling 
ofD^ = dN/duo^ the number of modes per unit uo^ is: 
Doj ~ uj^~^ in d dimensions. However, this scaling as- 
sumes that cj^ ~ at low frequency. On one hand, it 
has been pointed out that the anomalous low frequency 
dispersion of u'^ ^ q for low will change the expected 
scaling to ~ uj^ [25^ . On the other hand, experimen- 
tal data [7] and some hard-sphere Monte Carlo simulation 
results [23 for (111) seem to suggest that ~ u'^ at 
finite but low cj, in accord with the standard Debye pre- 
diction for a 3D material. We now compute the DOS, of 
finite rhombic (111) and square (100) patches cut from 
the full FCC crystal, paying special attention to the low 
frequency scaling behavior. 

We first compute the set of eigenvalues of ^, Ap where 
p indexes the mode. The p-th. frequency is then given 
by ujp = 1/ y/mXp where m is particle mass. We report 
frequency in units of ujq = ^/K/m. We first consider the 
integrated DOS: N{uj) = D^duj. Figure [g] shows N{uj) 
for (111) and (100) patches with P = 32 atoms along each 
edge embedded in a FCC crystal with C = 32 cubic unit 
cells along each edge. In addtion, we have overplotted the 
conventional Debye scaling uj ~ A/"^/^ and the expected 
modified Debye scaling: u ~ A/'^/^, for an intermediate. 
For both the patches, we see that the conventional Debye 
scaling does worse than the new modified-Debye predic- 
tion. However, N{uj) of the (100) patch shows almost 
perfect agreement with the expected non-Debye scaling 



(100) slice 
(111) slice 

t1/4 




10 



N 



100 



1000 



FIG. 6: Integrated DOS N{uj) = D^jduj. Frequencies, cj, 
as a function of index, N ^ for (100) and (111) FCC patches 
(P = 32, C = 32). 
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FIG. 7: Computed DOS, D^j^ normalized by the total num- 
ber of normal modes, N ^ for various FCC crystal patches. 
Top: Patch from (100) crystal plane. Bottom: Patch from 
(111) crystal plane. Upper panel: Doj/N scaled by the non- 
Debye scaling uo'^ . Lower panel: D^j/N scaled by the Debye 
scaling u;^. 



from A/" = 10 to about N = 200, whereas the (111) patch 
deviates by as much as 20% from the expected non-Debye 
scaling even in this low-cj regime. We verify that our re- 
sults are qualitatively unchanged if patch size P is kept 
constant while the embedding FCC crystal's size C is 
increased. That is, the size of the actual system has min- 
imal impact on the results, while the size of the obser- 
vation window must induce some finite size effects. We 
suspect that increasing the patch size would extend the 
modified-Debye scaling regime to lower frequencies, but 
that one would always observe departures by at least 20% 
for uo > O.GcjQ. 

We next make a histogram with the computed frequen- 
cies to directly estimate for various FCC patches and 
check scaling behavior at low uj. Figure [7| shows the DOS 
for (100) (top) and (111) (bottom) patches scaled by uj^ 
and o;^. In each case is normalized by the total num- 
ber of normal modes so that area under the graph is 
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unity. Our histogramming algorithm uses bins of vari- 
able width in cj-space so that the number of points lying 
in each bin is the same. For consistency we use 80 points 
per bin in all our histogram plots. We have checked that 
the location of the van Hove peaks and low frequency 
scaling behavior of is independent of these choices. 
We see from the plots of figure [7| that for C > P the 
embedding crystal is large enough such that C does not 
directly impact the results. The (100) patch DOS, when 
scaled by the expected non-Debye behavior, co'^, shows a 
plateau for low cj, while it lacks a plateau when scaled 
by the conventional Debye scaling, cj^; in accord with the 
integrated DOS in figure [6) 

In contrast to the (100) patch, the (111) behavior (fig- 
ure [t] bottom) is less clear. For the largest patches stud- 
ied, P = 64, there is a clear bump in the plateau in 
D^/uj^ at the lowest frequencies. This bump gives the 
conventionally scaled version, D^/uj'^, something close to 
an apparent plateau at low frequency; something that 
was clearly not observed in the (100) patch. An exper- 
imenter, using the statistical approach described later, 
looking at a patch of this size would be hard pressed to 
decide whether the low frequency regime scaled like the 
expected oo^ or the conventional o;^. 

In parting, we note that a similar plateau in 
D^/omega^ for a (111) patch was observed in refer- 
ence [23 . However, the DOS in that case was computed 
from a displacement correlation analysis in a hard sphere 
simulation. As we will show below, statistical artifacts 
themselves can act to give an even more robust plateau 
in / omego? at intermediate frequencies, so it becomes 
difficult to say whether in the case of reference [23 , the 
apparent plateau was the result of projection artifacts or 
statistical artifacts. 



C. Mode Structure 

We next focus our attention on the planewave decom- 
position of normal modes of the patch Green's function, 
and their relationship to the dynamic structure factor 
(DSF) [26], the average square amplitude of plane waves 
in thermal equiliubrium. We start by computing the full 
set of normal modes ^p, and eigenvalues \p of Q using the 
MATLAB eig() function, p indexes the modes. We then 
perform a spatial Fourier transform of each normal mode 
and split into longitudinal and transverse components, 
respectively, as follows: 



Next we compute the longitudinal (transverse) DSF, 
Sl{t)-, which is written in terms of normal modes as [27l 
28 : 



r 



(2) 



Here r labels a lattice site; E^{q) and E^{q) are square 
amplitudes of the longitudinal and transverse Fourier 
transform components, respectively, for a reciprocal lat- 
tice vector (of the patch) q. 



U — Urj 



muj 

p 

(3) 

Here m is particle mass and the temperature is such 
that hT/Ka^ = 1. The sum runs over modes in- 
dexed by p with corresponding energy eigenfrequency 
uOp = 1/^/raXp. Q {x) is the Heaviside step function, 
which is equal to 1 when x > and otherwise. S{q^ uS) , 
encodes the participation of various planewaves in a nor- 
mal mode of eigenfrequency uj [211. 

Figures and ^ show S[q^ uS) plotted in the first BZ for 
(111) and (100) FCC patches for a set of 3 values of uo that 
are typical low, mid-range and high eigenfrequencies of Q. 
We chose these frequencies so that the 'high' frequency 
is roughly near the second van Hove peak in the DOS, 
while the 'mid' frequency is near the first van Hove peak 
and the 'low' frequency is lower in the DOS. In both plots 
the upper row shows the longitudinal component of the 
DSF while bottom row shows the transverse component. 

The (111) patch DSF plots of figure [s] show that low 
energy normal modes (left column) have dominant par- 
ticipation from planewaves of a single wavenumber. The 
longitudinal planewaves participate at roughly half the 
wavenumber as transverse planewaves in modes of a given 
frequency uo and the DSF is isotropic as we expect from 
the low k dispersion of a (111) patch. For modes of in- 
termediate frequency (middle column) which correspond 
to the first van Hove peak in the DOS (see bottom 
plot of figure [t]), the dominant transverse planewaves 
are near the Brillouin Zone boundary. The DSF for 
higher frequency modes (right column) has peaks at 
multiple wavenumbers, however participating planewaves 
are still localized in the BZ. The dominant longitudi- 
nal planewaves are near the BZ boundary, as one would 
expect since uj is near the second van Hove peak. In 
addition, while the DSF has hexagonal symmetery, the 
longitudinal component is fairly isotropic except for the 
highest frequency normal modes. This is consistent with 
the pesudo-dispersion of Q of (111) patches seen in fig- 
ure [3l 

In contrast, the DSF for the (100) patch (fig- 
ure [9| shows participation from planewaves of different 
wavenumber in a given narrow frequency range. Roughly 
speaking, at a given cj, the DSF is peaked at the equi-cj 
contours of the dispersion plot. Much like the disper- 
sion of (100) patches, the DSF is highly anisotropic. At 
frequencies below the first van Hove singularity in the 
DOS (top plot of figure [7| the participating planewaves 
in the DSF are of relatively low wavenumber (left plot), 
whereas even for modes with frequency just beyond the 
first van Hove singularity we see dominant participat- 
ing transverse and longitudinal planewaves near the BZ 
boundary (middle and right column). 

We also compute an isotropic average of S{q^uj) in 
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FIG. 8: S[q,u) for a (111) FCC patch (P = 32, C = 32) for various u. Upper Row: Longitudinal. Lower Row: Transverse. 
Left: (jj/(jJo = 1.21. Center: uj/ujo = 1.65. Right: uj/ujo = 2.08. 
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FIG. 9: S[q,(jj) for a (100) FCC patch (P = 32, C = 32) for various cj. Upper Row: Longitudinal. Lower Row: Transverse. 
Left: cj/cjo = 1.21. Center: cj/cjq = 1.65. Right: cj/cjq = 2.08. 
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FIG. 10: Isotropically averaged DSF S{q,u), for a (111) 
FCC patch (P = 32, C = 32). The colorbar is log scale. Left: 
Longitudinal. Right: Transverse. 
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FIG. 11: Isotropically averaged DSF S{q,LS) for a (100) 
FCC patch (P = 32, C = 32). The colorbar is log scale. Top: 
Longitudinal. Bottom: Transverse. 
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and Inl show 6'(g,cj) computed for the 
Green's function oia (111) and a (100) patch, respec- 
tively. In each plot, the vertical axis spans the range of 
frequencies of Q while the horizontal axis spans the range 
of wavenumbers in the first BZ of the patch. 

In figure [lO] we see that for fixed, lower values of cj, 
S{q^uj) has a peak at a single value of q. Low energy 
modes of ^ for (111) patches have dominant participation 
by a planewaves of a single q. In contrast, high energy 
modes have multiple contributing planewaves. For exam- 
ple, modes around uj/ujq = 2.08, Sl{q,(^) (upper plot, 
right column of figure [8| have peaks for at least 2 values 
of q. As we expect from the pseudo-dispersion in fig- 
ure [Sj a longitudinally polarized planewave of wavenum- 
ber q participates in a normal mode of higher uj than a 



transverse planewave. In contrast to the (111) case, the 
isotropically averaged (100) DSF contains a wide range 
of q values at a given uj. This can be understood from 
the highly anisotropic dispersion. 



III. MONTE-CARLO SIMULATION 

Experimental studies have access to only finite statis- 
tical samples and will be subject to artifacts arising from 
incomplete statistical information. To study these ar- 
tifacts, we perform Monte-Carlo (MC) simulations using 
precisely the same Hamiltonian as our analytical calcula- 
tions. In particular we focus on the convergence of DOS. 

Using the Metropolis algorithm, we perform MC simu- 
lations to sample configurations from a Gibbs-Boltzmann 
distribution. Our system is a cubic FCC crystal 
with periodic boundary conditions, having C = 32 
cubic unit cells along each edge (4 x 32^ particles). 
The crystal Hamiltonian is constrained to be harmonic 
by specifying a Hessian Hiajf3 {i^j label atoms while 
label Cartesian axes) that includes only nearest- 
neighbor interactions. The energy is quadratic in par- 
ticle displacements {xia}, and is given by: E{{xia}) = 
^ ^iaj(3 ^ioiHiaj^Xjfs. We work in the same units as be- 
fore: energy is measured in units of Ko? . The tempera- 
ture is set so that k^T/Ka'^ = 1. We use MC steps whose 
size is sampled randomly from a uniform distribution of 
width 0.95a, where a is the nearest spacing. This gives 
us an acceptance ratio of approximately 81%. 

We assemble displacement covariance matrices to ob- 
tain the Green's function, Qiajis = {uiaUjfs) /k^T^ for 
(111) and (100) patches having P = 32 atoms along each 
edge. Here is the displacement of particle 'i' along 
axis a from its equilibrium position and the angle brack- 
ets denote averaging over MC sweeps (in one sweep we 
attempt to move each degree of freedom according to the 
Metropolis algorithm). We report results as a function 
of sampling time V, that we define as the number of 
(post equilibration) MC sweeps per degree of freedom in 
the patch (there are 2P^ degrees of freedom in a P x P 
patch) [16]. Note that r counts successive MC sweeps 
and does not address the issue of statistical independence 
between sampled configurations. 

We estimate the 'decorrelation time' in our MC setup 
by fitting the energy and single-particle displacement 
auto-correlation functions to an exponential function to 
obtain the characteristic decay time [29 . The decorrela- 
tion time is insensitive to size of the FCC crystal C. We 
find that the decorrelation time estimated from energy 
and displacement auto-correlation is roughly 8 and 24 
MC sweeps, respectively. These decorrelation times give 
some idea about the number of uncorrelated samples. 

In figure [12] we present the inferred frequency spectrum 
of Q of our (111) patch for various r together with the 
analytical result. We note that if r < 1 (the number 
of samples is less than the number of degrees of free- 
dom) we will necessarily pick up trivial modes (as we see 
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FIG. 12: Evolution of the frequency spectrum of the Green's 
function of a (111) FCC patch (P = 32, C = 32) inferred from 
Monte-Carlo simulations, r is the number of MC sweeps per 
patch degree of freedom. Top: Small sample. Bottom: Longer 
sample. 



FIG. 13: Scaled DOS, D^, oi 8. FCC crystal patch (P = 32, 
C = 32) normalized by the number of normal modes N — 
2P^, for different MC sampling times r. Top: Patch from 
(100) FCC plane. Bottom: Patch from (111) FCC plane. Up- 
per panel: Doj/N scaled by the non-Debye scaling cj^. Lower 
panel: D^j/N scaled by the Debye scaling cj^. 



in the spectrum corresponding to r = 0.977). Different 
parts of the frequency spectrum converge at quite differ- 
ent rates: the higher end of the frequency spectrum takes 
particularly long to converge. In fact, the change in the 
inferred spectrum with increasing r is so slow that if we 



only consider data in the top plot of figure [12] we might 
be tempted to conclude that the spectrum is converged 
at about r = 17.58 since the relative change per unit 
r in the spectrum is quite small. Indeed, small relative 
changes in the inferred spectrum are essentially the only 
tool an experimenter has to declare convergence since an 
analytical result isn't available [7l[T6l|l8]. The bottom 
panel of figure 12 shows the inferred spectrum and the 
analytical result in the limit of very large r. The data 
make it clear that convergence to the true spectrum is 
quite slow. 

Figure [13] shows of for various r. Again, we have 
scaled D^^ by the non-Debye prediction cj^, as well as the 
Debye law uJ^ . The low frequency behavior of the (100) 



DOS (top) is not severely impacted by the poor statistics 
and even for relatively small data sets, the D^juj^ data 
shows a much clearer plateau than D^/uj^. One would 
have little case to doubt the predicted non-Debye behav- 
ior. Surprisingly, the (111) patch (bottom plot) shows an 
apparent uJ^ scaling regime at intermediate uj that tends 
to disappear only at very large r. As we have seen in the 
previous section, even the true underlying spectrum, at 
least for these patch sizes, shows very little evidence for 
any uj^ scaling regime. So it appears that the artifacts 
from the finite statistics act in the same way as the arti- 
facts in the true underlying spectrum to push the result 
further away from the cj^ scaling toward an apparent uP' 
scaling. 

To make one final point, we plot the unsealed version 

In the main plot, 
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of the (111) DOS in Figure [14 
we normalize frequencies by the maximum observed fre- 
quency, ujrn^ in the data set as one might do with exper 
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FIG. 14: DOS, i:>a;, of a (111) FCC crystal patches (P = 32, 
C — 32) normalized by A/" = 2P^ for different MC sampling 
times r. Inset: D^o/N as a function of frequency cj. Main 
plot: Doj/N and cj, non-dimensionalized by the maximum 
frequency ujm- 



imental data, while the inset normalizes the frequencies 
by the underlying frequency scale in our known Hamilto- 
nian, ujq. We see that the inferred D^^ as a whole, shifts 
to lower frequencies relative to the entire spectrum for 
shorter MC sampling times. This is essentially due to the 
fact that ujrn shifts downward with increasing observa- 
tion time. Interestingly, at least with these normalization 
conventions on the energy scale, the effect of spurious in- 
ferred correlations arising from insufficient sampling act 
in a way very much like true underlying disorder. For 
infinitely good statistics, or no underlying disorder, one 
has sharp van Hove peaks in the DOS. As the statisti- 
cal sample is degraded, or one introduces disorder into 
the underlying model, the van Hove peaks broaden and 
shift to lower frequencies [15 . In effect, one may observe 
an excess of low frequency modes, or Boson peak, in the 
inferred spectrum because of insufficient statistics rather 
than any true underlying disorder. 



IV. DISCUSSION AND SUMMARY 

The ultimate use of the analysis in this paper is to 
inform experimenters studying glassy or crystalline col- 
loidal suspensions using particle trajectory data obtained 
from optical microscopy techniques. As we have shown, 
there are two issues that much be carefully addressed 
in analysis of such data: the effects of observing a re- 
stricted 2D portion of the system and imperfect statisti- 
cal information about the displacement correlations. We 
first described an analytical approach to understanding 
projection artifacts for systems governed by a harmonic 
Hamiltonian and then analyzed MC simulations based 
on the exact same Hamiltonian to study the statistical 



artifacts. The main utility of the approach we described 
here, and the advantage over previous work, is that both 
kinds of artifacts - projection and statistical - are ad- 
dressed within a single framework. 

As a side benefit, the symmetry properties of the pro- 
jected Green's function arise as a natural consequence of 
the symmetries of the embedding FCC lattice and the 
use of central forces. We did not need to adjust any free 
parameters describing the symmetry in the patch, e.g. 
as in reference [23 . A fully atomistic model also gives 
a realistic description of modes all the way up to the 
BZ boundary, whereas other approaches rely on ad- hoc 
discretization of an elastic continuum. 

We have shown that for any direction in the Bril- 
louin zone, the Green's function of a 2D patch embedded 
in an FCC crystal has an anomalous non-linear disper- 
sion reation cj^ ~ inline with earlier, long- wavelength 
limit, analytical prediction |221 [23] . Despite the same 
long wavelength scaling, the dispersion relation is very 
anisotropic for (100) patches, whereas it is fairly isotropic 
in the (111) case, at least at long wavelengths. Moreover 
we find that anisotropy is much more pronounced for 
the transverse dispersion of (100) FCC patches. Notably 
for (100) patches we find that along certain directions in 
the BZ, the longitudinal and transverse branches of the 
dispersion are nearly degenerate all the way to the BZ 
boundary. 

The non-linear, long wavelength, dispersion relation 
leads one to expect a non-Debye scaling, ~ o;^ in 
the DOS of the FCC patches at low frequencies. Our 
calculations show the (100) DOS has clear D^^ ~ cj^ scal- 
ing at low frequencies, however the picture is somewhat 
ambiguous for (111) patches and does not show a clear 
Doj ~ cj^ regime. This is similar to the observations from 
molecular dynamics simulations of FCC hard sphere sys- 
tems in reference ]23 . From the calculated DOS, that 
we show in figure [7| we see that the deviation from cj^ 
scaling at low uj can arise even with perfect statistical 
information. However, we do expect that this deviation 
itself is a finite (patch) size effect and that a ~ uo^ 
scaling regime would emerge for large enough patches. 

Next we studied the character of normal modes in the 
patches. The DSF gave the contribution of planewaves to 
normal modes of a given frequency. In both the (111) and 
(100) cases, the DSF at a given frequency was essentially 
determined by equi-o; contours of the dispersion. The 
strong anisotropy of the (100) dispersion made for broad 
contributions to the isotropically averaged DSF at a given 
uj due to the shape of the dispersion contour. 

The issue of statistical convergence was more subtle. 
We showed that the DOS converges surprisingly slowly 
with increasingly good statistical estimates of the dis- 
placement covariance. In general, we found that poor 
statistics smears out the van Hove singularities in the 
DOS and shifts the peak toward lower frequencies, much 
in the same way as true disorder. For the (111) patches, 
these sampling artifacts can give an apparent low fre- 
quency form for the DOS that appears to agree quite well 
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with the conventional Debye result. However, the DOS 
becomes closer to the expected non-Debye form with bet- 
ter statistics. 

In the future, we would like to extend our study to 
harmonic FCC crystals with springs drawn from a distri- 
bution (as in refs [HI [15]). It will be important to under- 
stand, quantitatively, the competing effects of incomplete 
statistical information and actual underlying disorder on 
the inferred spectrum. Beyond this, more realistic mod- 
els for interacting colloidal particles - such as hard sphere 
systems - bring anharmonicity and non- Gaussian behav- 
ior. Understanding all these effects in quantitative ways 
is necessary before precise connections can be made - in 
experiments or simulations - between the underlying in- 
teractions and the observed displacement correlations in 
small low- dimensional observation windows. 



Appendix 

A. FOURIER SPACE CALCULATION OF 
GREEN'S FUNCTION 

The Green's function of a crystal G = H~^, where 
H is its Hessian matrix. For a crystal that has a unit 
cell with a basis and a pair interatomic potential, we de- 
note the components of the Hessian matrix by: Hio^^j^y^ 
here z, j label unit cells, /i, v label atoms in the basis and 
label Cartesian axes. We consider periodic FCC 
crystals in a cubic cell with C, 4-atom cubic unit cells 
along each edge. This gives N = unit cells and 
M = AN atoms in the crystal. Translational invariance 
in the crystal implies that the Hessian is a function of 
only the separation Rij = Ri — Rj between unit cells, 
where Ri denotes the position of the ith. unit cell, thus 
giving Hiafj^j^u = Ha,^py{Rij) The set of allowed Rij is 
precisely the set of position vectors, R^ of the unit 
cells; suppressing indices labeling pairs of unit cells, we 
write the components of the Hessian as: Hafi/3u{R)' This 
allows us to represent H, which is a 3M x 3M matrix, 
in terms of 12 x 12 matrices H(^), whose compo- 
nents are Hoc^j3y{R). Note that having a 4-atom unit 
cell implies that: H(^) = H^(— ^), due to translational 
invariance of the lattice. 



Following Ashcroft and Mermin [26], we take a 3D 
Fourier transform of the Hessian matrix H(^) to com- 
pute the dynamical matrices according to: 



D(fe) 



R 



Here /c is a vector in the cubic reciprocal lattice and the 
sum runs over the set of unit cell position vectors R. 
The tilde denotes Fourier transformed operators. The 
dynamical matrices D(/c) are 12 x 12 complex, Hermitian 
matrices. Due to the identity H(^) = H^(— ^), we also 
have: D(-^) = D^(^). 

Next, we numerically invert the T>(k) to get the Fourier 
space Green's function G(^): 

Since there are 3 trivial translational modes of the crys- 
tal, G(0) is computed by taking a Moore- Penrose pseu- 
doinverse via a singular value decomposition (SVD). Fi- 
nally an inverse Fourier transform of the G(/c) gives us 
the complete real space Green's function G(^): 



G{R) 



1 



k 



In our implementation, we perform the forward Fourier 
transform as an explicit sum over k as very few of the 
H(^) are non-zero due to short range interactions; e.g. in 
the case of only nearest neighbor interactions there are 
at most 3^ = 9 non-zero H(^). However, the inverse 
Fourier transform is performed via a fast Fourier trans- 
form routine. 



Acknowledgments 

We thank Michael Widom for useful discussions at var- 
ious stages during this work. This material is based upon 
work supported by the National Science Foundation un- 
der Award Number NSF-DMR- 1056564. 



[1] R N. Pusey and W. V. Megen, Nature 320, 340 (1986). 
[2] U. Gasser, E. R. Weeks, A. Schofield, P. N. Pusey, and 

D. A. Weitz, Science 292, 258 (2001). 
[3] V. W. A. de Villeneuve, R. P. A. Dullens, D. G. A. L. 

Aarts, E. Groeneveld, J. H. Scherff, W. K. Kegel, and 

H. N. W. Lekkerkerker, Science 309, 1231 (2005). 
[4] E. R. Weeks, J. C. Crocker, A. C. Levitt, A. Schofield, 

and D. A. Weitz, Science 287, 627 (2000). 
[5] A. M. Alsayed, M. F. Islam, J. Zhang, P. J. CoUings, and 

A. G. Yodh, Science 309, 1207 (2005). 



[6] C. Brito, O. Dauchot, G. Biroli, and J.-P. Bouchaud, Soft 

Matter 6, 3013 (2010). 
[7] D. Kaya, N. L. Green, C. E. Maloney, and M. F. Islam, 

Science 329, 656 (2010). 
[8] K. Chen, W. G. Ellenbroek, Z. Zhang, D. T. N. Chen, 

P. J. Yunker, S. Menkes, C. Brito, O. Dauchot, W. van 

Saarloos, A. J. Liu, et al., Phys. Rev. Lett. 105, 025501 

(2010). 

[9] A. Ghosh, V. K. Chikkadi, P. Schall, J. Kurchan, and 
D. Bonn, Phys. Rev. Lett. 104, 248305 (2010). 



12 



[10] K. Zahn, A. Wille, G. Maret, S. Sengupta, and 

P. Nielaba, Phys. Rev. Lett. 90, 155506 (2003). 
[11] H. H. von Griinberg, P. Keim, K. Zahn, and G. Maret, 

Phys. Rev. Lett. 93, 255703 (2004). 
[12] Z. Cheng, J. Zhu, W. B. Russel, and P. M. Chaikin, Phys. 

Rev. Lett. 85, 1460 (2000). 
[13] A. J. Kurd, N. A. Clark, R. C. Mockler, and W. J. 

O'Sunivan, Phys. Rev. A 26, 2869 (1982). 
[14] W. Schirmacher, G. Diezemann, and C. Ganter, Phys. 

Rev. Lett. 81, 136 (1998). 
[15] S. N. Taraskin, Y. L. Loh, G. Natarajan, and S. R. El- 

hott, Phys. Rev. Lett. 86, 1255 (2001). 
[16] S. Henkes, C. Brito, and O. Dauchot, Soft Matter 8, 6092 

(2012). 

[17] G. F. Mazenko, Equilibrium Statistical Mechanics (John 

Wiley and Sons, Inc., 2000). 
[18] A. Ghosh, R. Mari, V. Chikkadi, P. Schall, J. Kurchan, 

and D. Bonn, Soft Matter 6, 3082 (2010). 
[19] K. Chen, M. L. Manning, P. J. Yunker, W. G. Ellenbroek, 

Z. Zhang, A. J. Liu, and A. G. Yodh, Phys. Rev. Lett. 

107, 108301 (2011). 
[20] A. Ghosh, V. Chikkadi, P. Schall, and D. Bonn, Phys. 

Rev. Lett. 107, 188303 (2011). 



[21] N. L. Green, D. Kaya, C. E. Maloney, and M. F. Islam, 

Phys. Rev. E 83, 051404 (2011). 
[22] C. A. Lemarchand, A. C. Maggs, and M. Schindler, EPL 

(Europhysics Letters) 97, 48007 (2012). 
[23] M. Schindler and A. Maggs, The European Physical 

Journal E: Soft Matter and Biological Physics 34, 1 

(2011) , ISSN 1292-8941, 10.1140/epje/i2011-11115-7. 
[24] M. Schindler and A. C. Maggs, Soft Matter 8, 3864 

(2012) . 

[25] A. Ghosh, R. Mari, V. Chikkadi, P. Schah, A. Maggs, 
and D. Bonn, Physica A: Statistical Mechanics and its 
Applications 390, 3061 (2011), ISSN 0378-4371. 

[26] N. Ashcroft and N. Mermin, Solid State Physics (Saun- 
ders College, Philadelphia, 1976). 

[27] H. Schober, Journal of Physics: Condensed Matter 16 
(2004). 

[28] H. Shintani and H. Tanaka, Nat Mater 7, 870 (2008). 
[29] M. Newman and G. Barkema, Monte Carlo Methods in 

Statistical Physics (Oxford University Press Inc., New 

York, 1999). 
[30] (100) results are qualitatively the same 



